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Abstract 

In  this  paper,  we  present  a new  technique  for  generating  error  equidistributing  meshes 
that  satisfy  both  local  quasi-uniformity  and  a preset  minimal  mesh  spacing.  This  is  firstly 
done  in  the  one-dimensional  case  by  extending  the  Kautsky  and  Nichols  method  [6]  and 
then  in  the  two-dimensional  case  by  generalizing  the  tensor  product  methods  to  altern- 
ating curved  line  equidistributions.  With  the  new  meshing  approach,  we  have  achieved 
better  accuracy  in  approximation  using  interpolatory  radial  basis  functions  (RBFs).  Fur- 
thermore improved  accuracy  in  numerical  results  have  been  obtained  for  a class  of  linear 
and  non-homogeneous  PDEs  solved  by  the  dual  reciprocity  method  (DRM). 


1 Introduction 

The  adaptive  mesh  algorithms  have  been  widely  used  in  the  numerical  solution  of  par- 
tial differential  equations  (PDEs)  for  boundary  value  problems  [1,  13].  One  undesirable 
feature  of  an  error  equidistributing  mesh  is  that  there  is  no  guarantee  of  it  being  suffi- 
ciently smooth.  For  our  applications  of  interpolation  (using  RBFs),  the  distance  between 
points  becoming  too  small  can  imply  that  the  underlying  interpolation  matrix  becomes 
ill-conditioned. 

In  this  paper,  we  propose  a method  to  deal  with  this  problem  in  Section  2.  Essentially 
our  method  consists  of  modifying  the  error  monitor  function  in  a suitable  way  and 
then  equidistributing  the  new  function  so  that  the  minimal  mesh  size  constraint  can  be 
satisfied.  We  deal  with  the  extension  of  adaptive  mesh  to  two  dimensions  in  Section  3. 
Finally,  some  numerical  results  will  be  given  in  Section  4. 

2 An  adaptive  mesh  with  minimal  mesh  size  control 

In  the  ID  case,  a typical  adaptive  mesh  problem  can  be  stated  as  follows:  given  a mesh 
(uniform  or  non-uniform)  to,ti,...,tm,  and  its  corresponding  error  values  (usually  es- 
timated from  the  numerical  solution  using  a monitor  function  [5])  /o,  /i, . . .,  fm,  we  wish 
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to  find  a new  mesh 


n . xq  , , . . . , xn , 


(2.1) 


that  is  locally  bounded  with  respect  to  a positive  constant  k > 1 such  that  l/k  < 
hj/hj-i  < k,-j  = 1,2, 1,  hj  = Xj+ 1 — Xj,  while  the  errors  are  equidistributed  on 
mesh  II.  One  solution  to  this  problem  was  given  in  [6]  by  replacing  fj  by  fj  followed  by  a 
standard  equidistribution  algorithm,  fj  is  referred  to  as  the  padded  function  and  the  main 
idea  of  replacing  fj  is  increasing  the  values  of  the  function  /,  where  too  small,  to  prevent 
considerably  large  mesh  sizes.  We  now  propose  a method  of  further  modifying  fj  in  such 
a way  that  the  resulting  equidistribution  mesh  satisfies  the  preset  minimal  mesh  size 
hmin ■ Before  proceeding,  we  consider  replacing  the  piecewise  linear  function  /( x)  (with 
endpoint  values  fj  = f(tj))  by  another  piecewise  linear  function  Z(x)  (with  endpoint 
values  Zj  = f(xj)).  This  is  a technical  approximation  to  simplify  the  presentation; 
actually  the  proposed  method  may  work  without  this  step.  Note  that  if  we  were  to 
equidistribute  Z(x),  the  resulting  mesh  would  not  differ  from  Xj  much;  define  the  average 
value  of  the  monitor  function  as 

dl  = d\Z)  = ^(Zj+Zj+1)^.  (2.2) 

3=  0 

Our  aim  now  is  to  modify  some  Zj  values  so  that  the  modified  average  value  is  the  same 
as  d!  while  the  modified  values  ensure  a preset  minimal  mesh  size  /imin  is  satisfied.  To 
present  our  method,  we  note  that  insisting  on  hj  > hmin  implies  Zj  < Z where 

Zhmin  — d (2*3) 


and  Z is  the  critical  constant  to  realize  hmin.  This  points  a way  of  modifying  those  large 
values  of  Zj.  However  it  is  not  obvious  how  to  ensure  the  new  and  modified  average 
values  are  the  same,  i.e.  equidistribution  is  maintained  for  the  same  error  constant. 
Suppose  that  among  the  current  Zj  values,  there  are  M + 1 of  them  that  are  larger  than 
Z (i.e.  whose  corresponding  mesh  size  is  less  than  hmin)',  denote  these  values  by  Zkj  for 
J = 0,1,...,  M.  This  means  that  Zkj  < Z for  j = M + 1,  M + 2, . . . , n.  Here  the  sequence 
ko,  ki, ..  .,kn  represents  a permutation  of  0, 1, 2, ..  ,,n. 

It  turns  out  that  a suitable  modification  (from  Zj  to  Zj)  is  the  following: 


(i)  Zkj  = Z when  Zkj  > Z,  i.e.  for  j = 0, 1, . . .,M, 

M 


(ii)  Zkj=zkj+ —^L. 


E ^ u 

l=M+ 1 


E(^  - z^hk> 


i= 0 


hki 


(2.4) 


for  j = M+  l,M  + 2,  ...,n, 


where 


hki  — 


{hki  + hki- 1)/2  when  ki  0 ,n, 

ho/2  when  ki  = 0, 

hn- 1/2  when  ki  = n. 


(2.5) 
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For  a simple  illustration,  see  the  plot  of  Fig  3b.  To  prove  that  the  above  modification  is 
suitable,  we  first  present  the  following  result  for  a simple  case. 

Theorem  2.1  Let  xo,  xi, . • -,xn  be  a non-uniform  mesh  with  the  mesh  sizes  h:)  = 
Xj+i  — Xj  and  Z$,  Z\, . . .,  Zn  are  the  corresponding  error  values.  If  the  critical  constant 
value  Z as  in  (2.3),  and  only  one  value  Z\>  Z (i.e.  M — 1 and  all  others  Zj  are  less 
than  or  equal  to  Z),  the  modification  (2.4)  takes  the  following  form, 

j (i)  Z0  = Z0,  Z,  = Z, 

| (ii)  Zj  = Zs  + [(Z,  - Z)(h0  + ft,)/ 2]  /(ft,  + hj- ,)/2  for  j = 2, 3 

Then  the  average  value  d = d(Z)  of  the  modified  values  Zj  is  the  same  as  d'  = d'(Z)  in 
(2.2). 

Note  M = 1 here;  in  fact  the  results  holds  for  any  one  value  Zj  > Z.  Now  we  are  ready 
to  present  the  main  result  on  equation  (2.4)  with  regard  to  minimal  mesh  size  control. 

Theorem  2.2  With  the  error  function  modified  as  in  (2.4),  the  new  mesh  hj  resulting 
from  equidistribution  satisfies  (i)  the  average  error  value  remains  as  d' ; (ii)  hj  > hmin. 
Here  hmin  cannot  be  specified  to  be  larger  than  h = 1 fn  (the  uniform  mesh  size); 
practically  we  found  hm,n  £ [h?,h/2]  is  adequate.  Full  proofs  to  these  results  will  be 
given  in  the  full  version  of  this  paper  [10]. 

In  the  method  in  (2.4),  the  values  of  Zk]  which  are  less  than  but  close  to  Z may 
become  unnecessarily  larger  (e.g.  larger  than  Z)  and  therefore  we  can  propose  a further 
refinement.  We  can  keep  some  of  the  Zkj  values  which  are  between  Zj 2 and  Z.  In  other 
words,  we  only  modify  the  very  large  and  very  small  values  of  Z/,j  (see  plot  of  Fig  3b). 
Then  our  theorems  are  still  valid  but  the  proofs  may  need  minor  changes.  Finally  we 
summarise  our  adaptive  method  with  minimal  mesh  size  control  as  follows  (see  the  plot 
of  Fig  3b  for  an  illustration). 

Algorithm  2.3.  (Numerical  algorithm)  For  given  non-uniform  mesh  a = to,  t\, . . 
tm,  = b,  the  error  values  /o>  /i , • • ■ , fm,  values  c and  hmjn: 

(1)  Does  the  locally  bounded  mesh  algorithm  converge  to  the  new  mesh  a = Xq  < X\  < 

• • • < xn  = b which  is  sub-equidistributing  with  respect  to  c and  f,  that  is,  for  a 
sufficiently  large  value  of  the  integer  n such  that  J/6  / < nc,  and  the  inequalities 

rxi+ 1 

/ / < c,  j = 0, 1 , ....  77  - 1 

Jxj 

are  satisfied. 

(2)  Check  the  minimal  mesh  size  and  compare  it  with  the  . If  it  is  less  than  hmin, 
go  to  the  Step  3 otherwise  stop. 

(3)  Approximate  the  padding  values  Zj  = f(xj)  corresponding  to  the  new  mesh  by  using 
piecewise  linear  interpolation  of  fj  values  and  calculate  the  average  value 

j 71  1 I 

~^2(zj  + zj+i)-f,  where  hj=xj+1-xj, 
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and  Z according  to  Zhmin  = d. 

(4)  Obtain  the  decreasing  arrangement  of  Zj,  Zkj  by  ordering  them. 

(5)  Modify  the  Z kj  values  as  follows, 

(i)  Zkj  = Z when  Zkj  > Z. 

assuming,  that  for  j = 0, Zkj  > Z, 

(ii)  Zkj  = Zkj  when  Z/2  < Zkj  < Z , 


assuming,  that  for  j 
zk. 


(Hi)  Zk.  — Zk.  y 

J J Z-,I=AT+1 

/or  j 


zk, 


j = N + 1, N + 2, . . ,,n, 


ZtL0(zki-z)hki  /hkj, 


where  hki  was  introduced  in  (2.5). 

(6)  Check  the  modified  values  Zkj  in  the  stage  (Hi)  of  the  Step  5.  If  Zkj  < Z / 2 for  all 
j,  go  to  Step  7 otherwise  repeat  Step  5. 

(7)  Perform  the  equidistribution  procedure  for  the  modified  values  Zkj  arid  obtain  the 
new  adapting  mesh. 


3 Extension  to  two  dimensions 

The  concept  of  adapting  mesh  in  one  dimension  is  well  known  (see  e.g.  [5,  3]).  Extension 
of  this  idea  to  two  dimensions  is  not  straightforward.  For  a given  function  f(x,y ) and 
2D  domain  Cl,  an  obvious  extension  is  dividing  the  domain  Cl  into  some  subdomains 
in  such  a way  that 


Fig.  1 . In  Fig  (a)  the  monitor  values  corresponding  to  the  new  mesh  are  represented  by 
**’,  the  linear  interpolation  for  these  values  is  shown  by  and  in  Fig  (b)  the  modified 
values  of  the  padded  function,  represented  by  dash  line,  are  compared  with  the  original 
values. 
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FlG.  2.  In  Fig  (a)  equidistribution  of  slabs  in  the  two  coordinate  direction  and  in  Fig 
(b)  three  stages  of  the  new  method  are  shown. 

But,  such  a partition  is  not  unique  and  furthermore  satisfying  condition  (3.1)  properly 
is  not  simple.  Consequently,  this  condition  has  to  be  replaced.  Among  the  methods 
given  to  satisfy  the  condition  (3.1)  as  much  as  possible,  two  well  known  methods  are 
transformation  and  dimension  reduction.  Transformation  methods  are  based  on  mapping 
the  physical  domain  into  a simple  domain  with  a uniform  mesh  and  ultimately  applying 
the  equidistribution  condition  to  obtain  an  adapting  mesh  in  the  physical  domain  [4,  12]. 
These  methods  are  generally  costly  and  complicated  in  theory.  In  this  work  we  first 
consider  the  latter  method  which  is  easier  and  cheaper  than  the  former  method.  We 
then  present  a new  technique  to  generate  a 2D  mesh. 

3.1  Dimensions  reduction 


We  assume  that  D is  a rectangle  in  the  form  ft  = {{x,y),  a < x < b,  c < y < d).  A 
simple  idea  is  to  produce  the  mesh, 


such  that 


a = x o < xi  < . ..<  x„_i  < xn  = b, 
c = 2/0  < 2/1  < • • • < 2/m— 1 < 2/m  = d , 


fx(x,  y)  dy  dx  = constant, 


(3.2) 


and 


fy(x,  y ) dx  dy  — constant, 


(3.3) 


where  fx(x,y ) and  fy(x,y ) are  the  monitors  in  the  x and  y directions  respectively  (see 
Fig  3.1a).  Obviously  the  generated  mesh  by  this  method  is  much  different  from  an 
equi-distributing  mesh  that  one  expects  from  (3.1).  Another  method  which  leads  to  a 
non-rectangular  grid  is  dimensional  splitting  [11].  We  now  describe  a new  method  of 
type  dimension  reduction. 
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Fig.  3.  In  Fig  (a)  the  mesh  generated  by  the  new  method  for  function  in  (3.6)  and  in 
Fig  (b)  the  resulting  mesh  when  restricting  the  minimal  mesh  size  as  hmin  = hj 2 for 
the  same  function  are  shown. 


3.2  A new  approach  for  a 2D  mesh 

The  idea  is  based  on  the  tensor  product  method  and  therefore  a non-rect angular  grid. 
We  start  with  a uniform  mesh  in  a rectangular  region  Q and  perform  the  method  in 
three  stages.  In  the  first  stage,  the  error  equidistributing  is  performed  for  each  line  in 
the  horizontal  direction  (see  the  first  part  of  Fig  3.1b),  that  is, 
rxi+ 1 

fx{x,  Vi)  dx  = constant  for  i = 0, 1, . . . , m.  (3.4) 

Jij 

In  the  next  stage,  the  mesh  is  redistributed  in  the  vertical  direction  along  the  new  grid 
lines  (see  the  second  part  of  Fig  3.1b),  that  is, 

r«i+i 

fy  (xj  ,y)dy  = constant  for  j = 0, 1, . . . , n, 


r 

Jxi 


J 

J Si 


(3.5) 

where  Sj+i  — s*  is  the  distance  between  two  consecutive  points  ( Xj,yi ) and  (xj,yi+i) 
along  the  new  lines.  In  the  final  stage,  equidistributing  is  repeated  in  the  horizontal 
direction  along  the  grid  lines  (the  last  part  of  Fig  3.1b).  One  can  observe  that  repeating 
this  procedure  usually  leads  to  a convergent  mesh.  According  to  our  experiments,  the 
number  of  iterations  to  achieve  convergence  is  at  most  five.  The  resulting  mesh  by  this 
procedure  for  function 

u(x,y)  = e(4~x2~4y2)  (3.6) 

when  applying  the  arc-length  monitor  is  shown  in  Figure  3a.  The  idea  of  controlling 
the  mesh  size  can  also  be  applied  in  this  technique.  The  generated  mesh  for  the  same 
function  when  the  mesh  sizes  are  restricted  to  hmin  = ft/2,  where  h is  the  mesh  size  in 
the  case  of  uniform  mesh,  is  given  in  Figure  3b. 

4 Numerical  examples 

In  this  part  the  affect  of  adapting  the  mesh  on  the  accuracy  of  interpolation  and  the  DRM 
is  considered.  In  the  following  examples,  the  infinity  norm  has  been  used  to  measure  the 


Fig.  4.  The  resulting  mesh  when  using  the  new  method  for  function  in  Examples  1 and 
2 are  shown  in  Figures  (a)  and  (b)  respectively. 


Method 

stage 

Function  (El) 

Derivative 

Function  (E2) 

Derivative 

uniform  mesh 

— 

5.1E-2 

9.5E-1 

1.3E-2 

2.2E-1 

Adaptive  mesh 

first 

5.4E-3 

1.6E-1 

2.5E-3 

1.3E-2 

with  control 

second 

5.4E-3 

3.0E-1 

2.1E-3 

1.0E-1 

third 

3.8E-3 

3.0E-1 

3.7E-3 

1.0E-1 

Adaptive  mesh 

first 

1.4E-2 

9.9E-2 

2.5E-3 

1.5-2 

without  control 

second 

2.2E-2 

7.5E-1 

2.1E-3 

1.0E-1 

third 

1.8E-2 

6.0E-1 

4.5E-3 

1.2E-1 

Tab.  1.  The  interpolation  error  for  Examples  1 — 2 using  adaptive  mesh  with  and 
without  control  the  mesh  sizes. 


accuracy,  that  is,  if  u and  u are  the  exact  and  approximate  values  respectively  then  the 
error  is  calculated  as 

eu  = ! |xi(m)  - u(x)||oo  = max  | u{x)  - u(aj)|. 

x£D 

A polynomial  RBF,  1 + r3,  has  been  employed  in  this  work. 


Example  4.1  We  check  the  interpolation  in  terms  of  the  RBFs  for  the  function, 

u(x,  y)  = (l  — e3l~3)  sin(1.5  tt  y),  (4.1) 

in  a rectangular  domain.  The  generated  mesh  for  this  function  is  shown  in  Figure  4a. 

Table  4 shows  the  affect  of  adapting  mesh  on  the  interpolation  accuracy  with  and  without 
controlling  the  mesh  sizes.  As  one  can  observe,  using  the  adapting  mesh  considerably 
improves  the  accuracy  in  comparison  with  the  case  of  uniform  mesh.  Moreover,  the  result 
in  the  case  of  controlling  the  minimal  mesh  size  is  better. 
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Example  4.2  In  this  example  we  first  check  the  function  f2(x , y)  = 0.5-0. 5 tanh(-4+ 
16  x2  + 16  y2)  and  then  solve  the  linear  PDE:  XJ2u  + + x/j^  + xyu  = d,  with  the 

Dirichlet  boundary  condition  over  the  elliptic  domain  x 2 + 4 y2  = 4,  where  d is  a known 
function  such  that  the  exact  solution  is  u(x,y)  = f2(x,  y). 

Again  from  Table  4,  we  see  improved  approximation.  We  apply  the  DRM  method  [7]  for 
solution,  where  the  domain  integrals  are  approximated  by  using  RBF  interpolation.  The 
adaptive  mesh  for  this  function  is  given  in  Fig.  4b  and  has  been  observed  to  give  rise  to 
improved  DRM  solution. 

5 Conclusions 

We  considered  a new  algorithm  for  producing  a locally  bounded  mesh  with  a preset 
minimal  mesh  size.  Such  a mesh  is  used  to  overcome  the  ill-conditioning  problems  asso- 
ciated with  radial  basis  function  interpolation.  Extension  of  the  idea  to  the  2D  case  is 
also  considered.  Some  preliminary  and  improved  numerical  results  are  given. 
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